rm(list=ls())
# Set working directory here
# NOTE: Everything will be loaded from & saved to this directory
setwd("")

# Load programs
source("procedures.R")
library(ggplot2)

# Random seed
set.seed(222)
# Sample size
n <- 1e+6
# Parameters to be saved here
thetas <- rep(NA, 6)
# Choose DGP
dgp.nr <- 4

# Graphics formatting parameters
mccolors <- function(...) gray.colors(..., end = 0.7) 
sel <- c(1, 3, 5, 6)
cols <- rep(c("gray70", "gray40", "gray20"), each = 2)[sel]
ltys <- rep(c(1, 2, 3), 2)[sel]
lwds <- c(3, 3, 6, 3, 3, 3)[sel]
nms <- c("Log", "Brier", "Spherical", "Boosting", "As1", "As2")[sel]
colb <- c(1, cols)
ltyb <- c(1, ltys)
lwdb <- c(5, lwds)
nmb <- c("True", nms)

# Simulate data
dat <- sim(dgp = dgp.nr, n = n)

# Load objective functions
of <- getof()

# Estimate parameter under various scoring rules
for (j in 1:6) thetas[j] <- optimize(of[[j]], y = dat$y, x = as.matrix(dat$x), interval = c(-10, 10))$minimum

# Save parameters to file
write.table(thetas, paste0("thetas", dgp.nr, "_", n, ".csv"), col.names = FALSE, row.names = FALSE)

# Draw prediction plot
pdat <- cbind(dat$p, sapply(thetas, function(s) logit(dat$grid*s)))

pdf(paste0("pred", dgp.nr, "_", n, ".pdf"), family = "Bookman")
matplot(x = dat$grid, y = pdat[, c(1, sel + 1)], type = "l", col = colb, lty = ltyb, 
        xlab = "x", ylab = "", lwd = lwdb)
dev.off()

# Draw classification plot
pdf.options(family = "Bookman")
pdf(paste0("class", dgp.nr, "_", n, ".pdf"), family = "Bookman")
x <- 2
c <- as.matrix(seq(from=0.01, to=0.99, by=0.01))
tr <- (dat$p.fct(x) > c)
gd <- sapply(logit(x*thetas), function(s) s > c)
matplot(cbind(tr, gd[, sel]),type = "l", x=c, col=colb, ylab="", xlab=expression(c), 
        lty=ltyb, lwd = lwdb)
dev.off()

# First legend
pdf("legend1.pdf", family = "Bookman", width = 10, height = 3)
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend("center", inset = 0, nmb,
       col = colb, lty = ltyb, lwd = lwdb, bty = "n", horiz = TRUE, cex = 0.8)
dev.off()

# Probability of different prediction
pd <- data.frame(sr = character(), c = factor(), val = numeric(), stringsAsFactors = FALSE)
srlabs <- c("Log", "Brier", "Spherical", "Boosting", "As1", "As2")

for (i in 1){
  for (j in (i+1):6){
    for (c in c(0.6, 0.7, 0.8, 0.9)){
      sr1 <- i
      sr2 <- j
      p1 <- logit(dat$grid * thetas[sr1])
      p2 <- logit(dat$grid * thetas[sr2])
      xd <- dat$grid[which(sign(p1 - c) != sign(p2 - c))]
      if (length(xd) > 0) aux <- (dat$cdf.x(xd[length(xd)]) - dat$cdf.x(xd[1])) else aux <- 0
      pd <- rbind(pd, data.frame(sr = as.character(paste(srlabs[sr1], "-", srlabs[sr2])), c = as.character(c), val = aux))
    }
  }
}
theme_set(theme_gray(base_size = 15))
pl <- qplot(y = reorder(sr, -order(sr)), x = val, shape = c, geom = "point", xlab = "", ylab = "", size = I(3.5),
      data = pd, xlim = c(0, 0.25))
ggsave(pl, file = paste0("diff", dgp.nr, "_", n, ".pdf"), family = "Bookman")

# Draw plot for share of correct predictions 
cgrid <- seq(from = 0.01, to = 0.99, by = 0.01)
x <- dat$x
true <- sapply(x, function(z) dat$p.fct(z) > cgrid)
share.correct <- matrix(0, length(cgrid), 6)
for (jj in 1:6){
  aux <- sapply(x, function(z) logit(z*thetas[jj]) > cgrid)
  share.correct[, jj] <- rowMeans(true == aux)
}

pdf(paste0("sc", dgp.nr, "_", n, ".pdf"), family = "Bookman")
matplot(x = cgrid, y = share.correct[, sel], type = "l", lwd = lwds, col=cols, ylab="", 
        xlab=expression(c), lty=ltys)
dev.off()     

# Second legend
pdf("legend2.pdf", family = "Bookman", width = 10, height = 3)
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend("center", inset = 0, nms, col = cols, lty = ltys, lwd = lwds, bty = "n", horiz = TRUE, 
       cex = 0.8)
dev.off()
